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We introduce four basic two-dimensional (2D) plaquette configurations with onsite cubic nonlin- 
earities, which may be used as building blocks for 2D PT-symmetric lattices. For each configuration, 
we develop a dynamical model and examine its VT symmetry. The corresponding nonlinear modes 
are analyzed starting from the Hamiltonian limit, with zero value of the gain-loss coefficient, 7. 
Once the relevant waveforms have been identified (chiefly, in an analytical form), their stability is 
examined by means of linearization in the vicinity of stationary points. This reveals diverse and, 
occasionally, fairly complex bifurcations. The evolution of unstable modes is explored by means of 
pH ' direct simulations. In particular, stable localized modes are found in these systems, although the 

O j. majority of identified solutions is unstable. 

PACS numbers: 63.20.Pw, 05.45.Yv, 03.75.Lm, 03.65.Ca, 11.30.Er, 02.40.Xx, 02.20.Sv 
0" 1 I. INTRODUCTION 

r\J \ The theme of VT (parity-time) symmetric systems was initiated in the works of Bender and collaborators [l| as 
^ ■ an alternative to the standard quantum theory, where the Hamiltonian is postulated to be Hermitian. The principal 
conclusion of these works was that PT-invariant Hamiltonians, which are not necessarily Hermitian, may still give rise 
to completely real spectra, thus being appropriate for the description of physical settings. In terms of the Schrodinger- 
type Hamiltonians, which include the usual kinetic-energy operator and the potential term, V(x), the PT- invariance 
admits complex potentials, subject to constraint that V*(x) = V{— x). 

Recent developments in optics have resulted in an experimental realization of the originally theoretical concept of 
the PT-symmetric Hamiltonians, chiefly due to the work by Christodoulides and co-workers Q (see also [H). It has 
been demonstrated that the controllable imposition of symmetrically set and globally balanced gain and loss may 
render optical waveguiding arrays a fertile territory for the construction of PT-symmetric complex potentials. The 
first two such realizations made use of couplers composed of two waveguides with and without loss [J (so-called 
passive VT— couplers), or, in more "standard" form, a pair of coupled waveguides, one carrying gain and the other 
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one loss In fact, more general models of linearly coupled active (gain-carrying) and passive (lossy) intrinsically 
nonlinear waveguides, without imposing the condition of the gain-loss balance, were considered earlier, and stable 
solitons were found in them @, including exact solutions Q (see also a brief review in Ref. [H). Recently, an 
electronic analog of such settings has also been implemented H, ■ Configurations with a hidden VT symmetry 
have been identified also in fine-tuned parameter regions of microwave billiards [ill ]. Effects of the nonlinearity in 
a Gross-Pitaevski equation on the VT properties of a Bose-Einstein condensate have been analyzed in [l^ . The 
possibility to engineer PT-symmetric oligomers (coupled complexes of a few loss-and gain-carrying elements) [l3| . 
which may include nonlinearity, was an incentive to a broad array of additional studies on both the few-site systems and 
entire PT-symmetric lattices [14l42l| . More recently, nonlinear PT-symmetric systems, incorporating PT-balanced 
nonlinear terms, have drawn considerable interest too [22j|-[25j]. 

Most of the PT-invariant systems considered thus far have been one-dimensional (ID) in their nature, although 



the stability of solitons in 2D periodic PT-symmetric potentials has also been recently investigated [26[ . Actually, 2D 
arrays of optical waveguides can be readily built [I?} (the same is true about other quasi-discrete systems, including 
electrical ones), hence, a natural question is whether PT-symmetric oligomers (and ultimately lattices built of such 
building blocks) can be created in a 2D form. This work aims to make a basic step in this direction, by introducing 
fundamental 2D plaquettes consisting, typically, of four sites (in one case, it will be a five-site cross). These configu- 
rations, illustrated by Fig. [I] arc inspired by earlier works on 2D Hamiltonian lattices described by discrete nonlinear 
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FIG. 1: (Color online) The different fundamental plaquette configurations (i.e., two-dimensional oligomers) including the linear 
balanced gain and loss. Among these, (a), (c) and (d) are PT-symmetric, while (b) is not in the strict sense, but it is interesting 
too, as an implementation of alternating gain and loss nodes in the plaquette pattern. The nodes are labeled so as to connect 
the gain-loss profiles to the evolution of individual nodes in dynamical simulations. The sets are coded by chains of symbols, 
with +, — and corresponding, respectively, to the linear gain, loss, or absence of either effect at particular sites. 



Schrodinger equations [28|, where diverse classes of modes, including discrete solitary vortices [2!| [3(J, have been 
predicted and experimentally observed [3ll |32| . The plaquettes proposed herein should be straightforwardly accessi- 
ble with current experimental techniques in nonlinear optics, as a straightforward generalization of the coupler-based 
setting reported in Ref. @ ■ We start from the well-established Hamiltonian form of such plaquettes in the conserva- 
tive form, gradually turning on the gain- loss parameter (7), as the strength of the VT- invariant terms, to examine 
stationary states supported by the plaquettes, studying their stability against small perturbations and verifying the 
results through direct simulations. Actually, in this work we focus on those (quite diverse, although, obviously, not 
most generic) modes that can be found in an analytical form, while their stability is studied by means of numerical 
methods. The analytical calculations and the manifestations of interesting features, such as a potential persistence 
past the critical point of the linear PT symmetry, are enabled by the enhanced symmetry of the modes that we 
consider below. It is conceivable that additional asymmetric modes may exist too within these 2D configurations. 

Our principal motivation for studying the above systems stems from the fact that realizations of PT-symmetry e.g. 
within the realm of nonlinear optics will be inherently endowed with nonlinearity. Hence, it is only natural to inquire 
about the interplay of the above type of linear systems with the presence of nonlinear effects. In addition to this 
physical argument, there exists an intriguing mathematical one which concerns the existence, stability and dynamical 
fate of the nonlinear states in the presence of PT-symmetric perturbations. In particular, previous works [ia m, m, 
l33l |34| point to the direction that neither the existence, nor the stability of PT-symmetric nonlinear states mirrors 
that of their linear counterparts (or respects the phase transition of the latter generically) . The presentation of our 
results is structured as follows. Section [TT| contains a part of the analytical results, including a detailed analysis of 
the PT— symmetry properties of the nonlinear Schrodinger type model, as well as the spectral properties of the linear 
Hamiltonian subsystems. Section IIIII is devoted to the existence, stability and dynamics of stationary modes in the 
nonlinear systems. Beside analytical results, it contains a detailed presentation of the numerical findings. In section 
ITVl we summarize conclusions and discuss directions for future studies. 



II. THE SETUP AND SYMMETRY PROPERTIES 



A. General techniques 

The dynamics of the 2D plaquettes that we are going to consider is described by a multicomponent nonlinear 
Schrodinger equation (NLSE) 

iu = H L u + H NL (u)u (1) 

built over a transposition-symmetric linear N x N Hermitian matrix Hamiltonian Ht, = -ffj and an additional 
nonlinear N x N matrix operator, Hnl(u) = -ffjvx( u )- To understand the symmetry properties of this NLSE, we 
first analyze the associated linear problem 

ill = H L u, (2) 

and check then whether the symmetry is preserved by the nonlinear term, 7?tvl(u)u. The analysis can be built, in a 
part, on techniques developed for other nonlinear dynamical systems with symmetry preservation [35l - l42l |. 

For the present setups, the time reversal operation T can be defined as the combined action of a scalar-type complex 
conjugation T, T 2 = /, and the sign change of time, t — > —t, in full accordance with Wigner's original work which 
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introduced these concepts [43|. For the linear Schrodinger equation ([2j and its solutions 



N 

u(t) - ^e-^ 4 u„, 

n=l 

Hlu u = E n u n , 



(3) 



this implies 



T(*5tu) 
iftT(u) 



T(iJ L u), 
HiT(u), 



Tu(t) = Tu(t)| t 



n=l 



-iE n t- 



(4) 



where the overbar denotes complex conjugation. From the actual form of the gain-loss arrangements in the 2D 
plaquettes we can conjecture the existence of certain plaquette-dependent parity operators V, with V 2 = I, which will 
render the Hamiltonians VT— symmetric, [VT, Hl] = 0. To find an explicit representation of these parity operators 
V, we use the following ansatz, 



V G 



together with the pscudo-Hcrmiticity condition 



niVxAf 



Hi 



[T,V] = 



VH L V. 



(5) 



(6) 



The latter follows trivially from Eq. ([5]), [VT,Hl] — 0, and the transposition symmetry, Hl = H\. These parity 
operators V will be used to check whether the corresponding nonlinear terms Ha/x(u) satisfy the same VT— symmetry 
In contrast to linear setups with the VT— symmetry being either exact ([VT,Hl] = 0, VTu oc u) or spontaneously 
broken ([VT,Hl] = 0, VTu u), the nonlinear setups considered in the present paper allow for sectors of exact 
■PT-Symmctry ([VT, H NL (u)} = 0, VTu oc u) and of broken "PT-symmetry (VTu ok u =^> [VT, H NL (u)} ^ 0), 
as it is common for nonlinear PT-symmetric systems [22l|-[25l|. 



B. PT— symmetry properties of 2D plaquettes 



We start from the 2D plaquettc of 0+0- type depicted as configuration (a) in Fig. [T] This plaqucttc has only 
two (diagonally opposite) nodes carrying the gain and loss, while the other two nodes bear no such effects. The 
corresponding dynamical equations for the amplitudes at the four sites of this oligomer are 



iua 
iiiB 
iiic 
Hid 



-k(u B + u D ) - \u A \ 2 u A , 

-k(iiA + uc) - \ub\ 2 ub + i^UB, 

-k(v,B + ud) - \uc\ 2 uc, 

-k(v,A + uc) - \ud\ 2 ud - i^uo, 



(7) 



where 7 G R is the above-mentioned gain-loss coefficient, and k £ R is a real coupling constant. The nonlincarity 
coefficients are scaled to be 1 (we use time t as the evolution variable, although in the mathematically equivalent 
propagation equations for optical waveguides t has to be identified with the propagation distance, z). 
Denoting u := (u A , ub,uc,ud) T € C 4 , the matrices Hl and Hnl{u) in (JT]) take the form of 



-k -k 



H L = 



H NL (u) = - 




(8) 



(9) 
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To find the parity matrix P which renders the linear Hamiltonian PT— symmetric, [PT, Hi] = 0, we use the pscudo- 
Hermiticity condition ([5]) and notice that 



Hi — Hifi - 
Hi.o = -k(I 
His 



H 



L,l, 



(10) 
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L,0> 
t 



-H 



LA' 



Obviously, the following relations should hold: 



VHi.oV = H Lfi 
VHr X V = -His 



[P,H Lfi ) = 
{T,H LA }=0. 



(11) 



The first of these conditions together with P 
one of the three types, 

Vox :=I®(Tx, 



I, P ^ I reduces the possible form of the parity transformation to 



PxO '■= &st 



(12) 



where cr x ,y,z are the usual Pauli matrices. Taking into account that o x g z o x = —o~ z , the anti-commutativity condition 
in Eqs. (|11[) singles out the only possible parity matrix: 



V = PxO = 



/ 
/ 



( 1 
1 
10 

\o 1 



(13) 



for the linear transformation, i.e., the matrix interchanging A and C, as well as B and D. One then immediately 
checks that 



H NL (Vu)=TH NL {u)T = - 



(\uc\ 



V o 





|U£>| 









M 









M 2 / 



(14) 



Hence, in contrast to the linear component Hl, the nonlinear terms Hni(u) corresponding to Eqs. ([7]) are not 
Vl~— symmetric, in the usual matrix sense. Rather, the symmetry properties of the nonlinear terms have to be 
considered in the context of the nonlinear Schrodinger equation itself. Acting with PT on Eq. ([!]) we observe that 



PT(id t u) = PT[H L u + H NL (u)u\, 
id t (PTu) = H L PTu + H NL {PTu){PTu). 



(15) 



Hence, the full NLSE system ([7]) remains invariant if we define the PT transformation of the vectorial wave function 
obeying this system as follows: 



■PTU 

Vu(-t) 



J4> 



u(t). 



(16) 



This is in full analogy to the condition of exact PT— symmetry for the corresponding linear Schrodinger equation 1 . 
But the condition of spontaneously broken PT— symmetry ([PT,Hi] = 0, PTu u) is replaced by the condition of 
completely broken PT— symmetry (PTu qk. u => [PT, Hnl(vl)] ^ 0). In contrast to the present 2D plaquettes, 
which are mainly motivated by feasible experimental realizations, one can envision more sophisticated setups with 
VH NL {u)P = Hlf L (u). This will lead to a new type of partial (or intermediate) PT— symmetry (to be considered 
elsewhere), which for solutions u(i) with broken PT— symmetry will keep the nonlinear term Hnl(u) explicitly 



1 ) We note that apart from trivial stationary type solutions with factorizing structure u(t) = e~' Et uo, nonlinearity matrices Hnl(u) 
of more general type than that in J9} and (1141 may be envisioned which may produce PT— symmetric solutions u(t) with less simple time 
dependence. A detailed analysis of such systems will be presented elsewhere. 
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Vl~— symmetric {V— pscudo-Hcrmitian) in the matrix sense, but not VT— symmetric (under inclusion of the explicit 
time reversal t — > —t) in the sense of the NLSE system. 

For configurations (b) and (c) in Fig.Q] -£/l,o and Hnl(u) are still given by Eqs. ([SJ and (|T0|) . but with 
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(17) 



(18) 



respectively. Hence, relations pT|) arc valid for both configurations (b) and (c) as well. Using (fT2"|) in VH^^V = —Hl,i, 
we find a richer variety of parity operators V than for configuration (a). Configuration (b) allows for 



V Qx 
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whereas configuration (c) may be associated with 
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For configuration (d), we have 
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(23) 



and simple computer algebra gives again two possible parity operators: 



V d n = 



/ 1 \ 

' 1 1 

10 

1 

V o o o i o J 



/ 1 \ 

' 1 ' 

10 

10 

V i o o o o / 



(24) 



in strong structural analogy to configuration (b). One verifies that Hatl('PTu) = 'PHnl(u)'P ^ H^ NL (u) holds also 
for configurations (b), (c) and (d), hence all 2D plaquettes considered in the present paper are not VT— symmetric 
in the usual matrix sense. 



C. Spectral behavior of associated linear setups 



Next, we turn to the eigenvalue problems of the linear setups associated with plaquettes (a) - (d), i.e., to solutions 
of the equation 

H L u n = E n u n . (25) 
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From the corresponding characteristic polynomials, dct(7?i — EI) = 0, 



(a) 
(b) 
(c) 



0. 



we find 





(d): 


E{E 2 + 


(a): 


El,2 


= 0, 


(b): 


El, 2 


= ±n> 


(c): 


E h2 


= y/2k 2 - 




E3,4 


= -^2k- 


(d): 


E\,2 


= 



E 2 [E 2 - (4fc 2 - 7 2 )] = 0, 
(E 2 +1 2 )[E 2 -(4k 2 --y 2 )] 
E 4 -2{2k 2 -~f 2 )E 2 + 7 4 = 0, 

1 2 )[E 2 -{Ak 2 - 1 2 )]=Q, 



(26) 



±y/l~k 2 



E 3A = ±^IW 



j 2 ±2k v r k 2 



- 7 2 ± 2k^k 2 - 7 2 
E 3A = ±^4fc 2 -7 2 , E 5 = 0. 



(27) 



Obviously, the matrix Hamiltonians Hl for plaquettes (a) and (d) are not of full rank. For plaquette (a) we find 
rank (if/,) = 2, and Hl has a two-dimensional kernel space, ker(iJ^) = span c (ui, U2). For plaquette (d) we find 
rank (iiz,) = 4 and ker(ffi) = C* x 115, where C* = C — {0}. Moreover, we see that the spectrum for plaquette (d), 
up to the additional eigenvalue E5 = 0, coincides with that for (b). The different eigenvalues of the 4-node plaquettes 
displayed in Eqs. ([27)) show that these plaquettes are also physically not equivalent. Equivalence classes of nonlinear 
4-node plaquettes with isospectral linear Hamiltonians Hl but different pairwisc couplings have been considered, 
e.g., in [3j|. For plaquettes (a), (b) and (d) an exceptional point (EP) occurs at 7 2 = 4fc 2 , being associated with a 
branching of the eigenvalue pair E 3A 



£3.4 G I 
£3,4 G il 



for 
for 



Ak 2 > 



7 



4fc 2 < 7 2 . 



(28) 



In the case of plaquette (a), all four eigenvalues are involved in the branching at 7 = ±2fc, where E\ = . . . = E4 = 0. 
Via Jordan decomposition (e.g., with the help of the corresponding linear algebra tool of Mathematica) we find that 



(a): 



H L (l = ±2k) 



( 1 
10 


Vo 



= J 3 (o)e Ji(o), 



(29) 



i.e., a spectral degeneration of the type (O^O 1 ) in Arnold's notation 44|, or, in other words, a third-order EP with 
a single decoupled fourth mode. Hence, plaquettes of type (a) may serve as an easily implementable testground for 
the experimental investigation of third-order EPs (see e.g. |45l448j ). For plaquettes (b) and (d) we have second-order 
EPs at 7 2 = 4fc 2 , similar as for plaquette (c) where a pair of second-order EPs occurs at 7 2 = k 2 with E\ = E2 = \k\, 
E 3 =E i = -\k\. 

From the eigenvalues in ([27)) we read off the VT— symmetry content of the four types of plaquettes. The sector of 
exact VT— symmetry (i.e., the sector with all eigenvalues purely real, E n G R, Vn) corresponds to 



(a) 
(b) 
(c) 
(d) 



7 2 < 4fc 2 , 
7 = 0, 
7 2 < k 2 , 
7 = 0, 



(30) 



i.e., for plaquettes (b) and (d) the VT— symmetry is spontaneously broken as soon as the gain- loss coupling is switched 
on, namely for 7 7^ 0. 

III. EXISTENCE, STABILITY AND DYNAMICS OF NONLINEAR STATES 



In this section, we seek stationary solutions of the type 

Uo(t) = e- lEt u a , £gK, u = (a, b, c, d) T e C 4 



(31) 
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constructed over constant vectors Uo. According to (|15|) and (|T6|) such solutions will be 'PT— symmetric provided it 
holds VTuq = e tip UQ for some <p G K. We will test these symmetry properties for the solutions to be obtained. We 
note that restricting the explicit analysis to stationary solutions of the type (j3~Tj) we by construction exclude from this 
analysis "PT— violating solutions with E ^ R which are necessarily non-stationary. 

A useful technical tool to facilitate the explicit derivation of stationary solutions Uo(t) are conservation equations 
of the type 



d t (u^Yu) = mt (iffy - YH) u, 



(32) 



constructed from Eq. ([T]) and its adjoint, where Y denotes an arbitrary constant matrix. The most simplest of them 
can be found via Eqs. © and Eq. (fTUf to be 



5 t |u| 2 =5 t (ut u ) = -2mt^ iilU 



= iu) 



Hl L {u)V-VH NL {u) u 



(33) 
(34) 



For stationary equations uo(£) the time-dependent phase factors e~ lEt cancel so that the left-hand-sides of these rela- 
tions vanish, yielding simple algebraic constraints on the right-hand-sides. From Eq. (|34p we see that for stationary so- 
lutions the VT inner product 2 will remain conserved {xvVm =const) regardless of the violated V— pseudo-Hermiticity, 



T'H 1 NL (u)'P iJjvx(u), characteristic for of our specific nonlinear plaquette couplings (see Eq. (fl"4"|) ). 

Subsequently, we first derive classes of stationary solutions Uo(t) explicitly. Then, we analyze the stability of small 
perturbations over these stationary solutions by the linearization, via ansatz 



u(t) 



,-iEi 



u + S(e xt r + e M s) + 0(S 2 ) 



xi , 



\S\ « 1, 



(35) 



where 8 is the small amplitude of the perturbation. Exponents A can be defined as Wick-rotated eigenvalues from the 
corresponding 8x8 perturbation matrix B (see, e.g., [28[ for more details): 



(B-iA/ 8 )x = 

td . ( 9 Un F(u,u) d an F(u,u) \ „_i 9 q 4 
{-d Un F(u,u) -d Un F(u,u)J ' 

x = (r,5) T , (36) 

where 

F(u, u) := [H(u) - E]u (37) 

characterizes the stationary problem, and the elements of the matrix B are evaluated at u = uo. Linear stability is 
ensured for A G iK, whereas A ^ iM corresponds to growing and decaying modes, i.e., exponential instabilities. 



A. The plaquette of the 0+0- type 

Substituting ansatz (|3"Tj) for the stationary solutions in Eqs. (fTJ), (f3"3")) and P4]l we obtain the following algebraic 
equations: 



Ea = 


fc(6H 


-d) 


+ \a\ 2 a, 




Eb = 


k(a - 


he) 


+ \b\ 2 b- 


ijb, 


Ec = 


fc(feH 


-d) 


+ |c| 2 c, 




Ed = 


k(a - 


he) 


+ \d\ 2 d^ 


-ijd, 



2 ) For completeness, we note that in the context of VT quantum mechanics (PTQM) the VT inner product was introduced first by 
Znojil in [49ll in 2001. Immediately afterwards, it was interpreted by Japaridze as indefinite inner product [5tJ in a Krein space and 
generalized by Mostafazadeh to the r\— metric in the context of pseudo-Hermitian Hamiltonians l5lll . Finally, it was used by Bender, 
Brody and Jones in 2002 to construct the positive definite CVT inner product [52l [. For oligomer settings (of plaquettes or other few 
site configurations), it can be employed, e.g., to derive a simple algebraic constraint or as a criterion of the numerical accuracy of the 
evolutionary dynamics (especially since the solutions rapidly acquire very large amplitudes when unstable, as will be seen below). It also 
turned out useful in [33ll , 
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c> t |u| 2 = -2ivJH Ltl u, 

= 2 7 (|6| 2 -|d| 2 ), (39) 

and 



d t (vSVu) = iut H^ NL {u)V-VH NL (u) 



u. 



= (\a\ 2 -\c\ 2 )(ac-ca) + (\b\ 2 -\d\ 2 )(bd-db). (40) 
These equations can be analyzed via the Madclung substitution (i.e., via amplitude-phase decomposition), 

Ae^" , b = Be 14 " 1 , c = Ce 1 ^ , d = De 1 ^ . (41) 



a 



Without loss of generality, we may fix <p a = 0. 

For arbitrary phase factors in (|41|) . Eqs. (|39| and (|40[) are satisfied by A — C and B = D. Using this condition in 
Eq. p8[) and dividing each equation Q38p by the phase factor on its left-hand side, one obtains the imaginary parts of 
the resulting equations: 






= kB 


[sin( 


(fib - 


K)- 


h sin( 


(fid - 0a)] 


7B 


= kA 


sin( 


fia~ 


<M- 


- sin( 


fic - <f>b)} 





= kB 


[sin( 


(fib - 


h)- 


- sin( 


fid - 0c)] 


7 B 


= kA 


sin( 


fia - 


M- 


- sin( 


fic - 4>d)] 



(fib + (fid J_\ I (fib - (fid 

COS 



COS 



2 



(fib + (fid , \ ( (fib — (fid 
COS 



2 V 2 

0a + 0c 



d cos 



(42) 



For ^> a = the first of these equations implies sin(0 o ) = — sin(0<j), hence either 0j = — 0d (case 1) or (fid = (fib — 7r 
(case 2). In case 1, we conclude from the third equation that either (j) b ^ ±7r/2 and C = (case la), or <\> b = ±ir/2 
and 4> c is arbitrary (case lb). In case 2 the third equation is satisfied automatically. In all the three cases, the second 
and the fourth equation are compatible. They give 

case la: sin(0 6 ) = -^j-^, 0c = 0, <j) d = -(fib, 

-yB 

case lb: cos(0 c ) = =F— - 1, (fid = -(fib = Tn/2, 

■yB 

case 2: sin(0 6 ) + sin(0 fc - <j> c ) — —rr > (fid = (fib -n. (43) 

Returning to the phase-factor divided equations (|38p and considering their real parts, we find 

EA = kB [cos((j) b - c/) a ) +cos((t) d - c/) a )}+ A A , 
EB = kA[cos((t) a -(t) b ) + cos((t) c -(f> b )}+B 3 , 
EA = kB [cos((j> b - 4> c ) + cos(<p d - <f> c )] + A 3 , 

EB = fcA[cos(0 a -0 d )+cos(0 c -0 d )]+S 3 . (44) 
The pairwise compatibility of the first and third, as well as of the second and fourth equations requires 

COs(0 b - 0a) + COs((f> d - (f) a ) = COs((/) b - <f> c ) + COs((f> d - <j> c ), 

cos(0 a - (f> b ) + cos(0 c - <j) b ) = cos(0 a - (j) d ) + cos(0 c - <f> d ). (45) 
For case la, these conditions are trivially satisfied, whereas for the remaining cases they lead to further restrictions: 

2kA 



case lb: <j) c = 0; tt ==> 7 = T 



B 

7 = 0; 

jB 

case 2: 4> c = 2<f> b ± n, sin(0 b ) = -^T^- ( 46 ) 



9 



In this way the phase angles are fixed for all the three cases and we can turn to the amplitudes. The corresponding 
equation sets reduce to 

case la: EA = 2kB cos(4> b ) + A 3 , 

EB = 2kAcos(cj) b ) + B 3 , 
case lb,2: E = A 2 = B 2 . (47) 

In the latter two cases (lb and 2) the amplitudes and phases completely decouple and we have 

A = B = C = D = y/\E~\. (48) 

Case la allows for a richer behavior. Equating the terms 2kcos((f>t) in the upper two equations (|4"T|) leads to the 
constraint 

A 2 (E - A 2 ) = B 2 (E - B 2 ), (49) 
which can be resolved by A = B (case laa) as well as by E = A 2 + B 2 (case lab). The analysis of these two cases 

I j2 B 2 

can be completed with the help of the relation cos(0b) = ±\/ 1 ttTo from Eq. (1431) . 

V 4k z A z 

As result we obtain the following set of stationary solutions: 

7# 

case la: sin(0 fc ) = -^r-p ^ = 0, 4> d = -<f>b, 



case laa: A = B = C = D=\jE^ ^jAk 2 - 7 2 , (50) 

2kA 

case lab: A = C, B = D = ■ E = A 2 + B 2 , (51) 

V^ 4 + 7 2 

case lb: 4>d = — 4>b = T 7T /2, c = O,7r, 7 = ±2k, 7 = 0, 

A = B = C = D = VE, (52) 

7 

case 2: sin(0f,) = -— , 4> d — §b ~ tt, 4>c — ^4>b ± 7r, 



A = B = C = J D = v£. (53) 

From Eq. (|40|) . it can also be seen that either A = C or if A ^ C, then sin(0 a — (j> c ) = must be true. Here, we 
use the information available so far to check the VT— symmetry content of the solutions (|50]) - ([53| explicitly. For 
stationary solutions u(t) = e~ lEt uo, E <G K the PT— symmetry condition (fT6|) implies 

TTu = e^u . (54) 



Taking into account that T acts as complex conjugation, we see from the explicit structure of V = V x q in Eq. f| 13[) 
that a stationary solution is VT— symmetric if, with <j) a = 0, it has <j) c = and <j)d = —<f>b (up to a common phase 
shift). Additionally, the amplitudes have to coincide pairwise: A = C, B = D, For Eqs. ([50)) - (f53l this means 
that all case-1 stationary solutions with <j> c = are PT— symmetric in their present form. The case-2 mode becomes 
explicitly VT— symmetric after a global U(V) multiplication by a phase factor: 

u = e l (^-"/ 2) v , 

V := A \ e - i (4>t,-K/2)^ e iK/2^ e i(4, b -ir/2)^ e -in/2'] T ^ 

VTv = v , (55) 



where cf> c = 2</>& — tt has to be chosen in Eq. (f53[) . We note that this procedure is effectively equivalent to a redefinition 
of the original phase constraint: <fi a = n- <fi a = — 4> b + tt/2 at the very beginning of the calculations in Eq. (|4"Tj) . 

The linear stability analysis was performed numerically. Subsequently we present corresponding graphical results. 
The plaqucttcs (b) - (d) can be analyzed in a similar way. For brevity's sake, in Fig. IIII Al we present only the basic 
numerical results, by means of the following symbols: 



• case laa with A = B = C = D= yE + y 4fc 2 + 7 2 - blue circles 
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Y Y 

FIG. 2: (Color online) Profiles of the solutions for plaquette (a) from Fig. [TJ with E — 2 and k = 1. Four different branches 
of the solutions are denoted by blue circles, red crosses, black squares and green stars. The top left and right panel display, 
respectively, the squared absolute values of the amplitudes and phase differences between adjacent sites for the respective 
states. The bottom left and right panels show real (the instability growth rates) and imaginary (oscillation frequencies) parts 
of the eigenvalues produced by the linearization around the stationary states. The continuations are shown versus the gain-loss 
parameter 7. 



• case laa with A = B = C = D= \JE — yj Ak 2 ~ 7 2 — red crosses; 

• case lab — green stars; 

• case 2 — black squares; 

• Case lb is not depicted explicitly because it corresponds to point configurations without gain-loss (7 = 0) and 
to exceptional point configurations 7 = ±2fc. 

Figure 1111 Al presents the mode branches (their amplitudes, phases, and also their stability) over the gain-loss 
parameter 7, starting from the conservative system at 7 = 0. The same symbols are used in Fig. |3l which displays 
typical examples of the spectral plane (A r , A^) for stability eigenvalues A = A r + i\i of the linearization; recall that 
the modes are unstable if they give rise to X r 7^ 0. Explicitly we observe the following behavior. 



• case laa with A = B = C = D= \l E + \J 4fc 2 + 7 2 - blue circles 
According to Fig. [3j the present solution is stable. Notice that, although featuring a phase profile, it cannot 
be characterized as a vortex state (the same is true for some other configurations carrying phase structure). 
Interestingly, the relevant configuration is generically stable bearing two imaginary pairs of eigenvalues. 



• case laa with A = B = C = D= \JE — yj 4fc 2 — 7 2 — red crosses. 

Obviously, this kind of solutions as well as the previous one exist up to the exceptional point 7 = ±2fc of the 
"PT-symmetry breaking in the linear system, where the two branches collide and disappear (leave the stationary 
regime and become nonstationary) . As seen in Fig. [3l the present branch has two eigenvalue pairs which are 
purely imaginary for small 7, but become real (rendering the configuration unstable) at 7 = 1.49 and then 
7 = 1.73, respectively. Ultimately, these pairs of unstable eigenvalues collide at the origin of the spectral plane 
with those of the previous branch (blue circles). 

• case lab — green stars. 

This stationary solution has a number of interesting features. Firstly, it is the only one among the considered 
branches which has two unequal amplitudes. Secondly, it exists past the critical point 7 = ±2fc of the linear 
system, due to the effect of the nonlinearity (the extension of the existence region for nonlinear modes was 
earlier found in ID couplers [22[ and oligomers [l3l |33|). Furthermore, this branch has three non-zero pairs of 
stability eigenvalues, two of which form a quartet for small values of the gain-loss parameter, while the third 
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FIG. 3: (Color online) The stability plots for plaquette (a) from Fig. [T]with E = 2 and k = 1, for different values of 7. The 
notation for different branches is the same as in the previous figure. All branches are shown for 7 = 0.5, 7 = 1.2, 7 = 1.6, and 
7 = 1.9 (top left, top right, bottom left, and bottom right panels, respectively). 



is imaginary (i.e., the configuration is unstable due to the real parts of the eigenvalues within the quartet). At 
7 = 1.17, the eigenvalues of the complex quartet collapse into two imaginary pairs, rendering the configuration 
stable, in a narrow parametric interval. At 7 = 1.24, the former imaginary pair becomes real, destabilizing the 
state again, while subsequent bifurcations of imaginary pairs into real ones occur at 7 = 1.28 and 7 = 1.74 (at 
the latter point, all three non-zero pairs are real). Shortly thereafter, two of these pairs collide at 7 = 1.76 and 
rearrange into a complex quartet, which exists along with the real pair past that point. 

• case 2 — black squares. In contrast to all other branches, this one is always unstable. One of the two nonzero 
eigenvalue pairs is always real (while the other is always imaginary), as seen in Fig. [3j This branch also terminates 
at the exceptional point 7 = ±2fc, as relation sin (</>b) = —7/ (2k) cannot hold at I7I > \2k\. This branch collides 
with the two previous ones via a very degenerate bifurcation (that could be dubbed a "double saddle-center" 
bifurcation), which involves 3 branches instead of two as in the case of the generic saddle-center bifurcation, 
and two distinct eigenvalue pairs colliding at the origin of the spectral plane. 

By means of direct simulations, we have also examined the dynamics of the modes belonging to different branches 
in Fig. 21 The stable blue-circle branch demonstrates only oscillations under perturbations. This implies that, despite 
the presence of the gain-loss profile, none of the perturbation eigenmodes grows in this case. Nevertheless, the three 
other branches ultimately manifest their dynamical instability, which is observed through the growth of the amplitude 
at the gain-carrying site [B, in Fig. Ufa)] at the expense of the lossy site (D). That is, the amplitude of the solution 
at the site with the gain grows, while the amplitude of the solution at the dissipation site loses all of its initial power. 
Depending on the particular solution, passive sites (the ones without gain or loss, such as A and C) may be effectively 
driven by the gain (as in the case of the black-square-branch, where the site A is eventually amplified due to the 
growth of the amplitude at site B) or by the loss (red-cross and green-star branches, where, eventually, the amplitudes 
at both A and C sites lose all of their optical power). 
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(a) blue circle branch (b) red cross branch 




(c) black square branch (d) green star branch 

FIG. 4: (Color online) The perturbed evolution of different branches from Figs. MI Al and 151 at 7 = 1.9. Thin solid, thick 
solid, thin dashed, and thick dashed curves correspond to nodes A, B, C, D in Fig. QJa), respectively. In panel (b), the plots 
pertaining to sites A and C [see Fig. QJa)] overlap. Similarly, pairs of the plots for (A,B) and (C,D) overlap in (c), and for 
(A,C) they overlap in (d). 



B. The plaquette of the -| — | — type 

We now turn to the generalized (not exactly "PT-symmetric) configuration 3 featuring the alternation of the gain 
and loss along the plaquette in panel (b) of Fig.fU Indeed, the absence of 'PT-symmetry in this case is mirrored in the 
existence of imaginary eigenvalues in the linear problem of Eqs. (|27[) . as soon as 7 7^ 0. The corresponding nonlinear 
solutions (with E ^ R) are not covered by the stationary solution ansatz ([?!]) ■ Stationary solutions (with E € R) 
solely belong to dynamical regimes below the concrete VT— thresholds. Apart from the two VT— symmetry violating 
solutions, there should exist at least two stationary solutions which we construct in analogy to [cf. Eqs. (|38[) ] from 

Ea = k(b + d) + \a\ 2 a — 17a, 

Eb = k(a + c) + \b\ 2 b + ijb, 

Ec = k(b + d) + \c\ 2 c-ijc, 

Ed = k(a + c) + \d\ 2 d + ijd. (56) 

Substituting the Madelung representation (|41[) and setting A = B = C = D (for illustration purposes, we focus 
here only on this simplest case), we obtain 

sin((/>b - (f) a ) + sm(c/) d - <j> a ) 
cos(</> & - <j) a ) + cos((/)d - (j) a ) 



sin(</> 6 - (f>c) + sm(<j) d - 4> c ) 
cos(</) b - 4> c ) + cos(<p d - 4> c ) = 



2 
fc' 

E-A 2 



(57) 
(58) 



3 ) For the terminology concerning exact VT— symmetry, spontaneously broken VT— symmetry and completely broken VT— symmetry 
see the discussion of Eqs. II15H and 1161 . 
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Further, fixing 



0, Eqs. and Q55D yield 



sin 0t = sin i 



-2-, A 2 =E±JM? 
2k v 



7 2 - 



(59) 



Obviously, the solution terminates at point 7 = ±2k. Similar to what was done above, the continuation of this branch 
and typical examples of its linear stability are shown in Figs. [51 and HIIBl respectively. From here it is seen that the 
blue-circle branch, which has a complex quartet of eigenvalues, is always unstable. In fact, the gain-loss alternating 
configuration is generally found to be more prone to the instability. The red-cross branch is also unstable via a similar 
complex quartet of eigenvalues. This quartet, however, breaks into two real pairs for 7 > 1.5, and, eventually, the 
additional imaginary eigenvalue pair becomes real too at 7 > 1.74, making the solution highly unstable with three 
real eigenvalue pairs. The manifestation of the instability is shown in Fig. typically amounting to the growth of 
the amplitudes at one or more gain-carrying sites. 




FIG. 5: (Color online) The continuation of mode (|59p and its stability, supported by plaquette (b) in Fig. [T] for E — 2 and 
k = 1. 
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FIG. 6: Two typical stability plots for branch (|59[l . for E = 2, k = 1 and 7 = 1 and 1.8, respectively. 
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(a) blue circles branch 



(b) red crosses branch 



FIG. 7: The perturbed evolution of the modes of type (|59|l at 7 = 1 corresponding to the left panel of Fig. IIII Bl The plots 
pertaining to sites B and D [see Fig. QJb)] overlap in both panels. 



C. The plaquette of the -| — | type 

We now turn to the plaquette in Fig. QJc), which involves parallel rows of gain and loss. In this case, the stationary 
equations are 

Ea = k(b + d) + \a\ 2 a — i"fa, 
Eb = k(a + c) + \b\ 2 b- ijb, 
Ec = k(b + d) + \c\ 2 c + i-fc, 

Ed = k(a + c) + \d\ 2 d + ijd. (60) 

In this case too, we focus on symmetric states of the form oiA = B = C = D [see Eq. (|4"Tj) ]. which gives rise to two 
solutions displayed in Fig. IIII CI represented by the following analytical solutions: 



A 2 = E-k± ^k 2 - 



0, sin < 



sin <pd 



2 



(61) 
(62) 



A 2 = E + k±y/k 2 

4>a = 0, 4> b = 7T, 4> c 



7T, sm<p d 



2 

fe' 



(63) 
(64) 



The analysis demonstrates that the branch with the upper sign in Eq. (|6ip is always unstable (through two real pairs 
of eigenvalues), as shown by blue circles in Fig. IIII Cl On the other hand, the branch denoted by the red crosses, which 
corresponds to the lower sign in Eq. (|6T|) is stable up to 7 = 0.86, and then it gets unstable through a real eigenvalue 
pair. The black-squares branch with the upper sign in Eq. (|63|) is always stable, while the green-star branch with 
the lower sign in Eq. (|63[) is always unstable. At the lincar-'PT-symmctry breaking point 7 = k, we observe a strong 
degeneracy, since all the three pairs of eigenvalues for two of the branches (in the case of the blue circles, two real 
and one imaginary, and in the case of red crosses — one real and two imaginary) collapse at the origin of the spectral 
plane. On the other hand, the black-squares branch is always stable with three imaginary eigenvalue pairs, while the 
green-star branch has two imaginary and one real pair of eigenvalues. Between the latter two, there is again a collision 
of a pair at the origin at the critical condition, 7 = k. Direct simulations, presented for 7 = 0.5 in Fig. [51 demonstrate 
the stability of the lower-sign black-squares branch, while the instability of the waveform associated with the blue 
circles and the green stars leads to the growth and decay of the amplitudes at the sites carrying, respectively, the 
gain and loss. Notice that at the parameter values considered here, the red-cross branch is also dynamically stable as 
shown in the top right panel of Fig. |H1 



D. The plaquette of the -| — OH — type 

Lastly, motivated by the existence of known "cross" -shaped discrete- vortex modes in 2D conservative lattices, in 
addition to the fundamental discrete solitons [HI, l29j , we have also examined the five-site configuration, in which the 
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FIG. 8: (Color online) The characteristics of the mode of the +- 1 type, supported by plaquette (c) in Fig. [TJ and given 

in analytical form by Eqs. (|61|l - (|64|l . for E = 2 and k = 1. The blue circles correspond to the completely unstable branch 
with the upper sign in Eq. (|61[) , while the red crosses pertain to branch with the lower sign, which is stable at 7 < 0.86. 
The black-square and green-star branches correspond to the upper and lower sign in Eq. (|63[) . respectively. The former one is 
always stable, while the later one is always unstable. All four branches terminate at the critical point | — y | = \k\ of the linear 
■PT-symmetric system. 



central site does not carry any gain or loss, while the other four feature a "PT-balanced distribution of the gain and 
loss, as shown in panel (d) of Fig. [TJ Seeking for stationary states with propagation constant, G [instead of E in Eq. 
(jnij, as in this case we reserve label E for one of the sites of the 5-site plaquette in Fig. [JJd)], we get: 

Ga = kc + \a\ 2 a + ija, 

Gb = kc + \b\ 2 b - ijb, 

Gc = k(a + b + d + e) + \c\ 2 c, 

Gd = kc+ \d\ 2 d + i-fd, 

Ge = kc + \e\ 2 e — i^/e. (65) 

Similarly as before, we use the Madelung decomposition a = Ae l ' l>a ,b = i?e^ b ,c = Ce^^d = De*,e = Ed 1 ^' , cf. 
Eq. (HTJ, and focus on the simplest symmetric solutions with A = B = D = E. Without the loss of generality, we set 
4> c = 0, reducing the equations to 



C 2 {G-C 2 ) = 4A 2 (G-A 2 ), 
(kC) 2 = {-/A) 2 + (GA-A 3 ) 2 1 

4>a = -4>b = 4>d = —<f>e- (66) 

We report here numerical results for parameters G — 15, k = 1 (smaller G yields similar results but with fewer 
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FIG. 9: (Color online) The perturbed evolution of the four branches of the analytical solutions given by Eqs. H61[) -(|64 | ). which 
correspond to Fig. IHI CI with 7 = 0.5. 



solution branches). We have identified five different solutions in this case, see Figs. lIII Dl and HnDl for the representation 
of the continuation of the different branches, and for typical examples of their stability (the latter is shown for 7 = 0.1, 
0.5 and 0.95). There are two branches (green stars and black squares) that only exist at 7 < 0.13, colliding and 
terminating at that point. One of them has three real eigenvalue pairs and one imaginary pair, while the other branch 
has two real and two imaginary pairs. Two real pairs and one imaginary pair of green stars collide with two real pairs 
and one imaginary pair of black squares, respectively, while the final pairs of the two branches (one imaginary for 
the green stars and one real for the black squares) collide at the origin of the spectral plane. These collisions take 
place at 7 = 0.13, accounting for the saddle-center bifurcation at the point where those two branches terminate. On 
the other hand, there exist two more branches (red crosses and magenta diamonds in Fig. 1111 D[) . which collide at 
I7I = |fc|. One of these branches (the less unstable one, represented by magenta diamonds) bears only an instability 
induced by an eigenvalue quartet, while the highly unstable branch depicted by the red crosses has four real pairs 
(two of which collide on the real axis and become complex at 7 > 0.92). Last but not least, the blue circles branch 
does not terminate at 7 = ±fc, but continues to larger values of the gain-loss parameter, |-y | > |fc|. It is also unstable 
(as the one represented by the magenta diamonds) due to a complex quartet of eigenvalues. 

The dynamics of the solutions belonging to these branches is shown in Fig. [T3J For the branches depicted by black 
squares and green stars (recall that they disappear through the collision and the first saddle-center bifurcation at 
7 = 0.13), the perturbed evolution is fairly simple: the amplitudes grow at the gain-carrying sites and decay at the 
lossy ones, while the central passive site (C) stays almost at zero amplitude. For the other branches, the amplitudes 
also grow at the two gain-carrying sites and decay at the lossy elements, while the passive site may be drawn to either 
the growth or decay. 



IV. CONCLUSIONS AND FUTURE CHALLENGES 



In the present work, we have proposed generalizations of the one-dimensional 'PT-symmetric nonlinear oligomers 
into two-dimensional plaquettes, which may be subsequently used as fundamental building blocks for the construction 
of 'PT-symmetric two-dimensional lattices. In this context, we have introduced four basic types of plaquettes, three 
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FIG. 10: (Color online) The characteristics of the different branches of solutions in the case of the five-site plaquette (d) in 
Fig. [Tjare shown for G — 15 and k = 1. The branches represented by the chains of black squares and green stars terminate 
at 7 = 0.13. The branches depicted by red crosses and magenta diamonds terminate at 7 = 1 [i.e., at the exceptional point 
I7I = |fc|], while the branch formed by blue circles continues past that point. 
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FIG. 11: (Color online) Case examples of the spectral planes of the linear-stability eigenvalues for the different solution branches 
shown in the previous figure, for G = 15, k = 1, and 7 = 0.1 and 0.95 (from left to right). 



of which in the form of four-site squares. The final one was in the form of the five-site cross, motivated by earlier 
works on cross-shaped (alias rhombic or site-centered) vortex solitons in the discrete nonlinear Schrodingcr equation. 
Our analysis was restricted to modes which could be found in the analytical form, while their stability against small 
perturbations was analyzed by means of numerical methods. Even within the framework of this restriction, many 
effects have been found, starting from the existence of solution branches that terminate at the critical points of the 
respective linear 'PT-symmetric systems — e.g., in the settings corresponding to plaquettes (a) and (c) in Fig. [TJ The 
bifurcation responsible for the termination of the pair of branches may take a complex degenerate form [such as the 
one in the case of setting (a)]. Other branches were found too, that continue to exist, due to the nonlinearity, past 
the critical points of the underlying linear systems. In addition, we have identified cases [like the gain-loss alternating 
pattern (b) or the cross plaquette of type (d)] when the VT symmetry is broken immediately after the introduction 
of the gain-loss pattern. The spectral stability of the different configurations was examined. Most frequently, the 
stationary modes are unstable, although stable branches were found too [e.g., in settings (a) and (c)]. We have also 
studied the perturbed dynamics of the modes. The evolution of unstable ones typically leads to the growth of the 
amplitudes at the gain-carrying sites and decay at the lossy ones. It was interesting to observe that the passive sites, 
without gain or loss, might be tipped towards growth or decay, depending on the particular solution (and possibly on 
specific initial conditions). 

The next relevant step of the analysis may be to search for more sophisticated stationary modes (that plausibly 
cannot be found in an analytical form) , produced by the symmetry breaking of the simplest modes considered in this 
work, cf. Ref. (22j . The difference of such modes from the 'PT-symmetric ones considered in the present work is 
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FIG. 12: (Color online) The perturbed evolution for solutions belonging to different branches from Figs. UlTDl and Fig. IIIIDl 
at 7 = 0.1. In panel (a), the amplitudes at the different sites of plaquette (d) from Fig. [1] (A,B,C,D,E) are depicted as follows. 
A: the line around 10" 1 ; B: the right one of the two triangle-like (oscillating) curves; C: the line around 10 ; D: overlapped 
by A; E: the left one of the two triangle like curves. In panel (b), the amplitudes at sites A and D overlap with each other 
and correspond to the bottom curve which tends to 0, while the amplitudes at sites B, C, E eventually grow to a large value. 
Panels (c) and (d) represent the dynamical effect of the gain at sites B and E, and loss at sites A and D, while the curve for 
the amplitude at site C remains very close to zero, (e) A and D overlap with each other and correspond to the bottom curve, 
which tends to 0; B, C, E eventually grow to values ~ 40. B and D overlap with each other and C starts a little higher than 
those two. 



the fact that modes with the unbroken symmetry form a continuous family of solutions, with energy E depending on 
the solution's amplitude, sec Eq. (|3ip . This feature, which is generic to conservative nonlinear systems, is shared by 
■pT-symmetric ones, due to the "automatic" balance between the separated gain and loss. On the other hand, the 
breaking of the symmetry gives rise to the typical behavior of systems with competing, but not explicitly balanced, 
gain and loss, which generate a single or several attractors, i.e., isolated solutions with a single or several values of 
the energy, rather than a continuous family. A paradigmatic example of the difference between continuous families 
of solutions in conservative models and isolated attractors in their (weakly) dissipative counterparts is the transition 
from the continuous family of solitons in the usual NLSE to a pair of isolated soliton solutions, one of which is an 
attractor (and the other is an unstable solution playing the role of the separatrix between attraction basins, the stable 
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soliton and the stable zero solution) in the complex Ginzburg-Landau equation, produced by the addition of the 
cubic-quintic combination of small dissipation and gain terms to the NLSE [53| . As concerns the systems considered 
in the present work, in the context of the breaking of the PT symmetry it may also be relevant to introduce a more 
general nonlinearity, which includes "PT-balanced cubic gain and loss terms, in addition to their linear counterparts 
(cf. Refs. [23| and [22j]). Nevertheless, it should also be noted that the issue of potential existence of isolated solutions 
versus branches of solutions in "PT-symmetric systems is already starting to be addressed in the relevant literature 
(including in plaqucttc-typc configurations), as in the very recent work of [33| . 

Moreover, the present work may pave the way to further considerations of two-dimensional 'PT-symmetric lattice 
systems, and even three-dimensional ones. In this context, the natural generalization is to construct periodic two- 
dimensional lattices of the building blocks presented here, and to identify counterparts of the modes reported here in 
the infinite lattices, along with new modes which may exist in that case. On the other hand, in the three-dimensional 
realm, the first step that needs to be completed would consist of the examination of a PT-symmetric cube composed 
of eight sites, and the nonlinear modes that it can support. This, in turn, may be a preamble towards constructing 
full three-dimensional PT-symmetric lattices. These topics are under present consideration and will be reported 
elsewhere. 



Acknowledgments 

UG thanks Holger Cartarius and Eva-Maria Graefe for useful discussions. PGK gratefully acknowledges support 
from the National Science Foundation under grant DMS-0806762 and CMMI-1000337, as well as from the Alexander 
von Humboldt Foundation and the Alexander S. Onassis Public Benefit Foundation. PGK and BAM also acknowledge 
support from the Binational Science Foundation under grant 2010239. 



References 



[1] C. M. Bender and S. Boettcher, Phys. Rev. Lett. 80, 5243 (1998); C. M. Bender, S. Boettcher and P. N. Meisinger, J. 
Math. Phys. 40, 2201 (1999). 

[2] Z. H. Musslimani, K. G. Makris, R. El-Ganainy and D. N. Christodoulides, Phys. Rev. Lett. 100, 030402 (2008); K. G. 

Makris, R. El-Ganainy, D. N. Christodoulides and Z. H. Musslimani, Phys. Rev. A 81, 063807 (2010). 
[3] S. Klaiman, U. Giinther, N. Moiseyev, Phys. Rev. Lett. 101, 080402 (2008). 

[4] A. Guo, G. J. Salamo, D. Duchesne, R. Morandotti, M. Volatier-Ravat, V. Aimez, G. A. Siviloglou and D. N. 

Christodoulides, Phys. Rev. Lett. 103, 093902 (2009). 
[5] C. E. Riiter, K. G. Makris, R. El-Ganainy, D. N. Christodoulides, M. Segev, D. Kip, Nature Phys. 6, 192 (2010). 
[6] B. A. Malomed and H. G. Winful, Phys. Rev. E 53, 5365 (1996); H. Sakaguchi and B. A. Malomed, Physica D 147, 273 

(2000); W. J. Firth and P. V. Paulau, Eur. Phys. J. D 59, 13 (2010); P. V. Paulau, D. Gomila, P. Colet, N. A. Loiko, N. N. 

Rosanov, T. Ackemann, and W. J. Firth, Opt. Express 18, 8859 (2010); A. Marini, D. V. Skryabin, and B. A. Malomed, 

ibid. 19, 6616 (2011); P. V. Paulau, D. Gomila, P. Colet, B. A. Malomed, and W. J. Firth, Phys. Rev. E 84, 036213 (2011). 
[7] J. Atai and B. A. Malomed, Phys. Rev. E 54, 4371 (1996). 
[8] B. A. Malomed, Chaos 17, 037117 (2007). 

[9] J. Schindler, A. Li, M. C. Zheng, F. M. Ellis and T. Kottos, Phys. Rev. A 84, 04010 1 (2011). 
[10] H. Ramezani, J. Schindler, F. M. Ellis, U. Giinther, and T. Kottos. larXiv:1205. 18471 

[11] S. Bittner, B. Dietz, U. Giinther, H. L. Harney, M. Miski-Oglu, A. Richter, and F. Schafer, Phys. Rev. Lett. 108, 024101 
(2012). 

[12] H. Cartarius and G. Wunner, arXiv: 1203. 1885 (to be published in the present special issue). 
[13] K. Li and P. G. Kevrekidis Phys. Rev. E 83, 066608 (2011). 

[14] H. Ramezani, T. Kottos, R. El-Ganainy and D. N. Christodoulides, Phys. Rev. A 82, 043803 (2010). 
[15] A. A. Sukhorukov, Z. Xu and Yu. S. Kivshar, Phys. Rev. A 82, 043818 (2010). 

[16] M. C. Zheng, D. N. Christodoulides, R. Fleischmann and T. Kottos, Phys. Rev. A 82, 010103(R) (2010). 
[17] E. M. Graefe, H. J. Korsch and A. E. Niederle, Phys. Rev. Lett. 101, 150408 (2008). 
[18] E. M. Graefe, H. J. Korsch and A. E. Niederle, Phys. Rev. A 82, 013629 (2010). 

[19] Z. Lin, H. Ramezani, T. Eichelkraut, T. Kottos, H. Cao and D. N. Christodoulides, Phys. Rev. Lett. 106, 213901 (2011). 

[20] S. V. Dmitriev, S. V. Suchkov, A. A. Sukhorukov, and Yu. S. Kivshar, Phys. Rev. A 84, 013833 (2011) 

[21] S. V. Suchkov, B.A. Malomed, S. V. Dmitriev and Yu. S. Kivshar, Phys. Rev. E 84, 046609 (2011); S. V. Suchkov, S. V. 

Dmitriev, B. A. Malomed, and Y. S. Kivshar, Phys. Rev. A 85, 033835 (2012). 
[22] A. E. Miroshnichenko, B.A. Malomed, and Yu. S. Kivshar Phys. Rev. A 84, 012123 (2011). 



20 



F. Kh. Abdullaev, Y. V. Kartashov, V. V. Konotop and D. A. Zezyulin, Phys. Rev. A 83, 041805 (2011) 
D. A. Zezyulin, Y. V. Kartashov, and V. V. Konotop, Europhys. Lett. 96, 64003 (2011). 
Y. He, X. Zhu, D. Mihalache, J. Liu, and Z. Chen, Phys. Rev. A 85, 013831 (2012). 
S. Nixon, L. Ge, and J. Yang, Phys. Rev. A 85, 023822 (2012). 

F. Lederer, G. I. Stegeman, D. N. Christodoulides, G. Assanto, M. Segev, and Y. Silberberg, Phys. Rep. 463, 1 (2008). 
P. G. Kevrekidis The discrete nonlinear Schrddinger equation: Mathematical Analysis, Numerical Computations and Phys- 
ical Perspectives, Springer- Verlag (Heidelberg, 2009). 

B. A. Malomed and P. G. Kevrekidis Phys. Rev. E 64, 026601 (2001). 

D. E. Pelinovsky, P. G. Kevrekidis, D. J. Frantzeskakis, Phys. D 212, 1 (2005). 

J. W. Fleischer, G. Bartal, O. Cohen, O. Manela, M. Segev, J. Hudock, and D. N. Christodoulides Phys. Rev. Lett. 92, 
123904 (2004). 

D. N. Neshev, T. J. Alexander, E. A. Ostrovskaya, Yu. S. Kivshar, H. Martin, I. Makasyuk, and Z. Chen, Phys. Rev. Lett. 
92, 123903 (2004). 

D. A. Zezyulin and V. V. Konotop, Phys. Rev. Lett. 108, 213906 (2012). 

E. M. Graefe, J. Phys. A (contribution to the present special issue). 
V. I. Arnold, Comm. Pure Appl. Math. 29, 557 (1976). 

J. Moser, Comm. Pure Appl. Math. 29 727 (1976). 

M. Golubitsky and I. Stewart, Arch. Rat. Mech. Anal. 87, 107 (1985). 

C. Elphick, E. Tirapegui, M. E. Brachet, P. Coullet, and G. Iooss, Physica D 29, 95 (1987). 
J. Montaldi, M. Roberts and I. Stewart, Nonlinearity 3, 695 (1990). 

G. M. Chechin and V. P. Sakhnenko, Physica D 117, 43 (1998). 

A. Ferrando, M. Zacares, P. Andrees, P. Fernandez de Cordoba and J. A. Monsoriu, Optics Express 13, 1073 (2005). 
M. Zacares, M. Arevalillo-Herraez and S. Abraham, Comp. Phys. Comm. 181, 35 (2010). 

E. Wigner, Group Theory and Its Application to Quantum Mechanics of Atomic Spectra, (Academic Press, 1959). 
V. I. Arnold, Russ. Math. Surv. 26, # 2, 29 (1972). 

U. Giinther and F. Stefani, Czech. J. Phys. 55, 1099-1106 (2005); |math-ph/0506021~j 
W. D. Heiss, J. Phys. A: Math. Theor. 41, 244010 (2008). 

E. M. Graefe, U. Giinther, H. J. Korsch, and A. E. Niederle, J. Phys. A: Math. Theor. 41, 255206 (2008). 
G. Demange and E. M. Graefe, J. Phys. A: Math. Theor. 45, 025303 (2012) 



M. Znojil, Rendic. Circ. Mat. Palermo, Ser. II, Suppl. 72 (2004), 211 - 218, |math-p h/0104012 
G. S. Japaridze, J. Phys. A: Math. Theor. 35, 1709-1 718 (2002), |quant-p h/0104077l~ 
A. Mostafazadeh, J. Math. Phys. 43, 205-214 (2002), |math-ph/010700T 



C. M. Bender, D. C. Brody, and H. F. Jones, Phys. Rev. Lett. 89, 270401 (2002), |quant-ph /0208076 
B. A. Malomed, Physica D 29, 155 (1987). 



